#define force_at_element_point_4(s1, s2, s3, s4, n)\
{\
  vscale(s1, elmv0, temp);\
  vscale(s2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(s3, elmv2, temp1);\
  vinc(temp1, temp);\
  vscale(s4, elmv3, temp1);\
  vinc(temp1, temp);\
  vdec(temp, elmi->vertices[n]->rforce);\
}

//Mingchen rforce = rforce - s1*displ1 - s2*displ2 - s3*displ3
//Mingchen stiffness matrix * displacement
#define force_at_element_point_3(s1, s2, s3, n)\
{\
  vscale(s1, elmv0, temp);\
  vscale(s2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(s3, elmv2, temp1);\
  vinc(temp1, temp);\
  vdec(temp, elmi->vertices[n]->rforce);\
}

#define force_at_integration_point_4(n1, n2, n3, n4)\
{\
  vscale(n1, elmv0, temp);\
  vscale(n2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(n3, elmv2, temp1);\
  vinc(temp1, temp);\
  vscale(n4, elmv3, temp1);\
  vinc(temp1, temp);\
  tempx = temp.x; tempy = temp.y; tempz = temp.z;\
  cu = cos(u); su = sin(u); cv = cos(v); sv = sin(v);\
  scucv = Sign(cu * cv); scusv = Sign(cu * sv); ssu = Sign(su);\
  pcue1 = (cu == 0.0) ? 0.0 : pow(fabs(cu), e1);\
  pcve2 = (cv == 0.0) ? 0.0 : pow(fabs(cv), e2);\
  pcue11 = (cu == 0.0) ? 0.0 : pow(fabs(cu), e1 - 1.0);\
  pcve21 = (cv == 0.0) ? 0.0 : pow(fabs(cv), e2 - 1.0);\
  psve2 = (sv == 0.0) ? 0.0 : pow(fabs(sv), e2);\
  psue1 = (su == 0.0) ? 0.0 : pow(fabs(su), e1);\
  psue11 = (su == 0.0) ? 0.0 : pow(fabs(su), e1 - 1.0);\
  psve21 = (sv == 0.0) ? 0.0 : pow(fabs(sv), e2 - 1.0);\
  fsc += tempx * a1 * scucv * pcue1 * pcve2 +\
          tempy * a2 * scusv * pcue1 * psve2 +\
          tempz * a3 * ssu * psue1;\
  fa1 += tempx * sc * scucv * pcue1 * pcve2;\
  fa2 += tempy * sc * scusv * pcue1 * psve2;\
  fa3 += tempz * sc * ssu * psue1;\
  fe1 += tempx * a1e1 * scucv * pcue11 * pcve2 +\
          tempy * a2e1 * scusv * pcue11 * psve2 +\
          tempz * a3e1 * ssu * psue11;\
  fe2 += tempx * a1e2 * scucv * pcue1 * pcve21 +\
          tempy * a2e2 * scusv * pcue1 * psve21;\
}

#define force_at_integration_point_3(n1, n2, n3)\
{\
  vscale(n1, elmv0, temp);\
  vscale(n2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(n3, elmv2, temp1);\
  vinc(temp1, temp);\
  tempx = temp.x; tempy = temp.y; tempz = temp.z;\
  cu = cos(u); su = sin(u); cv = cos(v); sv = sin(v);\
  scucv = Sign(cu * cv); scusv = Sign(cu * sv); ssu = Sign(su);\
  pcue1 = (cu == 0.0) ? 0.0 : pow(fabs(cu), e1);\
  pcve2 = (cv == 0.0) ? 0.0 : pow(fabs(cv), e2);\
  pcue11 = (cu == 0.0) ? 0.0 : pow(fabs(cu), e1 - 1.0);\
  pcve21 = (cv == 0.0) ? 0.0 : pow(fabs(cv), e2 - 1.0);\
  psve2 = (sv == 0.0) ? 0.0 : pow(fabs(sv), e2);\
  psue1 = (su == 0.0) ? 0.0 : pow(fabs(su), e1);\
  psue11 = (su == 0.0) ? 0.0 : pow(fabs(su), e1 - 1.0);\
  psve21 = (sv == 0.0) ? 0.0 : pow(fabs(sv), e2 - 1.0);\
  fsc += tempx * a1 * scucv * pcue1 * pcve2 +\
          tempy * a2 * scusv * pcue1 * psve2 +\
          tempz * a3 * ssu * psue1;\
  fa1 += tempx * sc * scucv * pcue1 * pcve2;\
  fa2 += tempy * sc * scusv * pcue1 * psve2;\
  fa3 += tempz * sc * ssu * psue1;\
  fe1 += tempx * a1e1 * scucv * pcue11 * pcve2 +\
          tempy * a2e1 * scusv * pcue11 * psve2 +\
          tempz * a3e1 * ssu * psue11;\
  fe2 += tempx * a1e2 * scucv * pcue1 * pcve21 +\
          tempy * a2e2 * scusv * pcue1 * psve21;\
}                  

#define force_at_integration_point_4_taper_bend(n1, n2, n3, n4)\
{\
  vscale(n1, elmv0, temp);\
  vscale(n2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(n3, elmv2, temp1);\
  vinc(temp1, temp);\
  vscale(n4, elmv3, temp1);\
  vinc(temp1, temp);\
  forcex = temp.x; forcey = temp.y; forcez = temp.z;\
  if(forcex != 0.0 || forcey != 0.0 || forcez != 0.0)\
   printf("fx = %f, fy = %f, fz = %f\n", forcex, forcey, forcez);\
  cu = cos(u); su = sin(u); cv = cos(v); sv = sin(v);\
  scucv = Sign(cu * cv); scusv = Sign(cu * sv); ssu = Sign(su);\
  if (cu < 0.0 || cu > 0.0) {\
    pcue1 = pow(fabs(cu), e1);\
    pcue11 = pow(fabs(cu), e1 - 1.0);\
  }\
  else {\
    pcue1 = 0.0;\
    pcue11 = 0.0;\
  }\
  if (cv < 0.0 || cv > 0.0) {\
    pcve2 = pow(fabs(cv), e2);\
    pcve21 = pow(fabs(cv), e2 - 1.0);\
  }\
  else {\
    pcve2 = 0.0;\
    pcve21 = 0.0;\
  }\
  if (su < 0.0 || su > 0.0) {\
    psue1 = pow(fabs(su), e1);\
    psue11 = pow(fabs(su), e1 - 1.0);\
  }\
  else {\
    psue1 = 0.0;\
    psue11 = 0.0;\
  }\
  kxsu = tx * ssu * psue1 + 1.0;\
  kysu = ty * ssu * psue1 + 1.0;\
  pcue1cve2 = scucv * pcue1 * pcve2;\
  pcue1sve2 = scusv * pcue1 * psve2;\
  spsue1 = ssu * psue1;\
  spsue11 = ssu * psue11;\
  pb13 = pi * b1 * b3;\
  pb123 = pi * b1 * b2 * b3;\
  comb = pi * b3 * (ssu * psue1 + b2 / sa3);\
  scomb = sin(comb); ccomb = cos(comb);\
  fsc4 += forcex * (kxsu * a1 * pcue1cve2 +\
                   pb123 * scomb / (sc * sa3)) +\
          forcey * kysu * a2 * pcue1sve2 +\
          forcez * a3 * ssu * psue1;\
  fa14 += forcex * kxsu * sc * pcue1cve2;\
  fa24 += forcey * kysu * sc * pcue1sve2;\
  fa34 += forcex * pb123 * scomb / (sa3 * a3) +\
          forcez * sc * spsue1;\
  fe14 += forcex * (tx * a1e1 * spsue11 * pcue1cve2 +\
                    kxsu * a1e1 * scucv * pcue11 * pcve2 -\
	            pb13 * e1 * spsue11 * scomb)+\
          forcey * (ty * a2e1 * spsue11 * pcue1sve2 +\
                    kysu * a2e1 * scusv * pcue11 * psve2) +\
          forcez * a3e1 * spsue11;\
  fe24 += forcex * kxsu * a1e2 * scucv * pcue1 * pcve21 +\
          forcey * kysu * a2e2 * scusv * pcue1 * psve21;\
  ftx4 += forcex * spsue1 * sc * a1 * pcue1cve2;\
  fty4 += forcey * spsue1 * sc * a2 * pcue1sve2;\
  fb14 += forcex * ccomb;\
  fb24 += -forcex * pb13 * scomb / sa3;\
  fb34 += -forcex * pi * b1 * (spsue1 + b2 / sa3) * scomb;\
}

#define force_at_integration_point_3_taper_bend(n1, n2, n3)\
{\
  vscale(n1, elmv0, temp);\
  vscale(n2, elmv1, temp1);\
  vinc(temp1, temp);\
  vscale(n3, elmv2, temp1);\
  vinc(temp1, temp);\
  forcex = temp.x; forcey = temp.y; forcez = temp.z;\
  cu = cos(u); su = sin(u); cv = cos(v); sv = sin(v);\
  scucv = Sign(cu * cv); scusv = Sign(cu * sv); ssu = Sign(su);\
  if (cu < 0.0 || cu > 0.0) {\
    pcue1 = pow(fabs(cu), e1);\
    pcue11 = pow(fabs(cu), e1 - 1.0);\
  }\
  else {\
    pcue1 = 0.0;\
    pcue11 = 0.0;\
  }\
  if (cv < 0.0 || cv > 0.0) {\
    pcve2 = pow(fabs(cv), e2);\
    pcve21 = pow(fabs(cv), e2 - 1.0);\
  }\
  else {\
    pcve2 = 0.0;\
    pcve21 = 0.0;\
  }\
  if (su < 0.0 || su > 0.0) {\
    psue1 = pow(fabs(su), e1);\
    psue11 = pow(fabs(su), e1 - 1.0);\
  }\
  else {\
    psue1 = 0.0;\
    psue11 = 0.0;\
  }\
  kxsu = tx * ssu * psue1 + 1.0;\
  kysu = ty * ssu * psue1 + 1.0;\
  pcue1cve2 = scucv * pcue1 * pcve2;\
  pcue1sve2 = scusv * pcue1 * psve2;\
  spsue1 = ssu * psue1;\
  spsue11 = ssu * psue11;\
  pb13 = pi * b1 * b3;\
  pb123 = pi * b1 * b2 * b3;\
  comb = pi * b3 * (ssu * psue1 + b2 / sa3);\
  scomb = sin(comb); ccomb = cos(comb);\
  fsc3 += forcex * (kxsu * a1 * pcue1cve2 +\
                   pb123 * scomb / (sc * sa3)) +\
          forcey * kysu * a2 * pcue1sve2 +\
          forcez * a3 * ssu * psue1;\
  fa13 += forcex * kxsu * sc * pcue1cve2;\
  fa23 += forcey * kysu * sc * pcue1sve2;\
  fa33 += forcex * pb123 * scomb / (sa3 * a3) +\
          forcez * sc * spsue1;\
  fe13 += forcex * (tx * a1e1 * spsue11 * pcue1cve2 +\
                    kxsu * a1e1 * scucv * pcue11 * pcve2 -\
	            pb13 * e1 * spsue11 * scomb)+\
          forcey * (ty * a2e1 * spsue11 * pcue1sve2 +\
                    kysu * a2e1 * scusv * pcue11 * psve2) +\
          forcez * a3e1 * spsue11;\
  fe23 += forcex * kxsu * a1e2 * scucv * pcue1 * pcve21 +\
          forcey * kysu * a2e2 * scusv * pcue1 * psve21;\
  ftx3 += forcex * spsue1 * sc * a1 * pcue1cve2;\
  fty3 += forcey * spsue1 * sc * a2 * pcue1sve2;\
  fb13 += forcex * ccomb;\
  fb23 += -forcex * pb13 * scomb / sa3;\
  fb33 += -forcex * pi * b1 * (spsue1 + b2 / sa3) * scomb;\
}                  

